% T1_dGemric: Estimates parameters and errorbars for T1 recovery in an IR sequence 
% Used signal model: S = A*[ 1-B*exp(-R1*TI) + exp(-R1*TR) ] 
% Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)

% After executing a file browser pops up, where the image to be mapped
% needs to be selected. Every inversion time should have it's own folder,
% and the images in each folder should have the same file name.

% Two ROIs need to be drawn: First a ROI with only noise and second the ROI
% where the mapping has to be performed

% There is a switch 'roibased' for performing voxel/roi based mappings at the top of
% this script.

% Output are dicom images of the R1map, R1errormap, T1map and T1errormap, saved
% in the directory of the first image of the serie. 
% IMPORTANT: The values of the R2 and R2 errormap are multiplied with 10000 to get into sensible gray values.
% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.

% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007.

clear all

% switch for voxel/roi/roimean type of estimation
roibased = 2; % 0=voxel, 1=roi all voxels, 2=roi mean
 
% choose if you want a T1 or a R1 map overlay image
overlaytype = 0; % 0=T1, 1=R1

[file,chosendir]=uigetfile('*.dcm');
cd(chosendir);
cd ..
studydir=pwd;
cd ..
outputdir = chosendir; %save maps in dir of chosen image

cd(studydir)
d=dir;
nfol=length(d);
fitdata(1).folder=char(d(1+2).name);

%read in the image that is picked by the user
chosenimage = dicomread(dicominfo(fullfile(chosendir,file)));

%get the dicom info
for i=1:nfol-2   
    fitdata(i).folder=char(d(i+2).name); 
    info = dicominfo([studydir,'\',fitdata(i).folder,'\',file]);
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime+10000;
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
    fitdata(i).TI=info.RepetitionTime;
end

%determine sigma of the noise
imshow(fitdata(2).image,[]);
text(35,120,'Select a region in the background','Color','yellow')
M=roipoly;
close;

%read in data to fitdata struct
for i=1:nfol-2
    D=fitdata(i).image.*M;
    sdyd(i)=std(nonzeros(D));
    maxsd=max(sdyd);
end

clear D;

%conversion factor from standard to rayleigh distribution
CRLBsigma=1.527*maxsd;

%draw masks
for i=1:nfol-2
    if i==1
       imshow(fitdata(2).image,[]);
       text(45,250,'Select the region of interest','Color','yellow')
%when ROI is to be drawn here:
       M=roipoly;
       close;
       
%when ROI is imported or defined
% load('C:\Users\Henk Smit\Desktop\MRData\110926 Jasper 10 dGe Volunteers\qmap_2_test_1_murid\Postprocessing\Mapping_Masks\_non_ARC_0001_fem_post.mat')
% load('C:\Users\Henk Smit\Desktop\MRData\110926 Jasper 10 dGe Volunteers\qmap_2_test_1_murid\Postprocessing\Mapping_Masks\TIR_2100_0010_fem_post.mat')

%        load('C:\Users\Henk Smit\Desktop\MRData\120320QmapKnees\001\Postprocessing\Mapping_Masks\TIR_2100_0010_fem_post.mat')        
%        M=logical(maskimage./maskimage);
%        M=zeros(256,256);
%        M(140:153,125:160)=1;
%        imshow(M)
%        M=maskimage;

    end
    fitdata(i).image=fitdata(i).image.*M;

end

for k=1:nfol-2
    TI(k)=double(fitdata(k).TI);
    TR(k)=double(fitdata(k).TR);
end

%initialize&allocate
TI=TI';
TR=TR';
yd=zeros(size(TI,1),size(TI,2));

if ~roibased
    [row,column]=find(fitdata(2).image);
    nrvoxels=size(row,1);
end;

dims=size(fitdata(1).image);
T1=zeros(dims(1),dims(2));
R1=T1;
A=T1;
B=T1;
STDCRT1=T1;
STDCRR1=T1;
STDCRB=T1;
STDCRA=T1;
 
%set optimization options
opt = optimset('fminunc');
opt = optimset(opt,'Diagnostics','off','DerivativeCheck','off','LargeScale','off','gradObj','on','Display','off','MaxIter',200,'Hessian','off','TolFun',1e-12,'Tolx',1e-10,'MaxFunEvals',200);
s = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,100],'Startpoint',[10000, 0, 0.001],'Maxiter',100,'Display','off','TolFun',10^-5,'TolX',10^-5);
% f = fittype('a-(b*exp(-x*c))','options',s);
fiteq = (['a*(1-b*exp(-x*c))']);
f = fittype(fiteq,'options',s);
fitrange = 1:size(TI);

%logistics for roibased all voxels
 if roibased == 1
     yd=double(fitdata(1).image);
     for k=2:nfol-2
            yd=[yd double(fitdata(k).image)];
     end
     yd=reshape(yd,size(yd,1)*size(yd,2),1);
     TI=repmat(TI,1,size(yd,1)/(nfol-2));
     TR=repmat(TR,1,size(yd,1)/(nfol-2));
     TI=TI';
     TR=TR';
     TI=reshape(TI,size(yd,1),1);
     TR=reshape(TR,size(yd,1),1);
     [nonzerox,nonzeroy]=find(yd);
     yd=yd(nonzerox,1);
     TI=TI(nonzerox,1);
     TR=TR(nonzerox,1);
     nrvoxels=1;
     
% logistics for roibased mean of voxels
 elseif roibased ==2
    yd=mean(double(nonzeros(fitdata(1).image)));
         for k=2:nfol-2
            yd=[yd mean(double(nonzeros(fitdata(k).image)))];
         end  
        yd=yd';
        nrvoxels = 1;      
 end

 % estimation loop
    % voxelbased fit
for i=1:nrvoxels
    disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels(1))]);
        if ~roibased
            if fitdata(1).image(row(i),column(i)) > 0;
                for k=1:nfol-2
                    yd(k)=double(fitdata(k).image(row(i),column(i)));
                end
                
                [LS,ML] = T1IRAbsComputePar(yd,TI,0,CRLBsigma,opt,s,f, TR, fitrange);
                [CR] = T1IRAbsCramerRao(ML,TI,CRLBsigma, TR);

                A(row(i),column(i))=ML(1,1);
                B(row(i),column(i)) = ML(2,1);
                T1(row(i),column(i)) = 1/ML(3,1);
                R1(row(i),column(i)) = ML(3,1);

                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                else
                    STDCRA(row(i),column(i)) = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB(row(i),column(i)) = sqrt(CR(2,2));
                else
                    STDCRB(row(i),column(i)) = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1(row(i),column(i)) = sqrt(CR(3,3));
                    STDCRT1(row(i),column(i)) = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1(row(i),column(i)) = 5;
                    STDCRT1(row(i),column(i)) = 10000;
                end
            else
                A(row(i),column(i))= 0;
                B(row(i),column(i)) = 0;
                T1(row(i),column(i)) = 0;
                R1(row(i),column(i)) = 0;

                STDCRR1(row(i),column(i)) = 0;
                STDCRT1(row(i),column(i)) = 0;
            end
        %roibased fit
        else 
                [LS,ML]=T1IRAbsComputePar(yd,TI,0,CRLBsigma,opt,s,f, TR, fitrange);
                [CR]=T1IRAbsCramerRao(ML,TI,CRLBsigma, TR);
                A = ML(1,1);
                B = ML(2,1);
                T1 = 1/ML(3,1);
                R1 = ML(3,1);
                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA = sqrt(CR(1,1));
                else
                    STDCRA = 10000;
                end
                if(CR(2,2)>0 && ~isnan(CR(2,2)))
                    STDCRB = sqrt(CR(2,2));
                else
                    STDCRB = 10000;
                end
                if(CR(3,3)>0 && ~isnan(CR(3,3)))
                    STDCRR1 = sqrt(CR(3,3));
                    STDCRT1 = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                else
                    STDCRR1 = 5;
                    STDCRT1 = 10000;
                end
        end
            
end

%output message



%write result in dcm files. R1 values*10000 to get in a sensible range
if ~roibased
    dicomwrite(int16(10000*R1),[outputdir,filesep,file(1:end-4),'_R1map','.dcm'],fitdata(1).info);
    dicomwrite(int16(T1),[outputdir,filesep,file(1:end-4),'_T1map','.dcm'],fitdata(1).info);
    %dicomwrite(int16(A),[maskdir,'\',maskfile,'Amap','.dcm'],fitdata(1).info);
    %dicomwrite(int16(B),[maskdir,'\',maskfile,'Bmap','.dcm'],fitdata(1).info);
    dicomwrite(int16(10000*STDCRR1),[outputdir,filesep,file(1:end-4),'_R1errormap','.dcm'],fitdata(1).info);
    dicomwrite(int16(STDCRT1),[outputdir,filesep,file(1:end-4),'_T1errormap','.dcm'],fitdata(1).info);
    
    disp(['Estimation finished. Used signal model: S = A*[ 1-B*exp(-R1*TI) + exp(-R1*TR) ] '  ])
    disp(['Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)'  ])
    disp(['Median A = ' num2str(median(nonzeros(A))) ', median B = ' num2str(median(nonzeros(B))) ', median T1 = ' num2str(median(nonzeros(T1)))  ])
    disp(['Resulting maps are stored in ' outputdir ])

    
    if ~overlaytype
        %T2 overlay image. The "[# #]" values are the range of the T2 map.
        figure('Name','Overlay with T1 map (milliseconds)', 'position', [800 600 500 500]);
        sc(T1,[0  2000], jet,sc(chosenimage, [0 max(max(chosenimage))],gray),T1==0);
        colorbar('south','XColor','white');
    else
        %R2 overlay image. The "[# #]" values are the range of the R2 map.
        figure('Name','Overlay with R1 map (R1*10000, 1/milliseconds)', 'position', [800 600 500 500]);
        sc(R1*10000,[0 100], jet,sc(chosenimage, [0 max(max(chosenimage))],gray),R1==0);
        colorbar('south','XColor','white');
    end

end

%display plots if roibased
if roibased
    disp(['Estimation finished. Used signal model: S = A*[ 1-B*exp(-R1*TI) + exp(-R1*TR) ] '  ])
    disp(['Parameter A = fully relaxed signal, B = inversion efficiency, R1 = 1/T1  (See McKenzie et al. 2006)'  ])
    scatter(TI,yd,'*k')
    hold on;
    ti=floor(0):1:round(max(TI));
    yy=abs(ML(1,1)*(1-ML(2,1)*exp(-ti*ML(3,1))+exp(-(TR(1)-TI(1)+ti)*ML(3,1))));
    plot(ti,yy,'--r');
    
    axis([0 max(ti) 0 1.4*max(yd)]);
    
%   option to plot the LS estimate too
%   yy2=abs(LS(1,1)*(1-LS(2,1)*exp(-ti*LS(3,1))+exp(-(TR(1)-TI(1)+ti)*LS(3,1))));
%   plot(ti,yy2)
    disp(['Model = Abs(A - B*exp(-R1*TI)'])
    disp(['Results: value +- standard deviation:'])
    disp(['R1 = ' num2str(ML(3,1)) ' +- ' num2str(sqrt(CR(3,3)))])
    disp(['T1 = ' num2str(1/ML(3,1)) ' +- ' num2str(sqrt(CR(3,3))/(ML(3,1)^2))])
    disp(['A = ' num2str(ML(1,1)) ' +- ' num2str(sqrt(CR(1,1)))])
    disp(['B = ' num2str(ML(2,1)) ' +- ' num2str(sqrt(CR(2,2)))])
end;

    
